本发明公开一种综合血管分布和视盘外观特性的视盘投影定位方法,其步骤为:(1)采用掩膜操作提取感兴趣视网膜眼底图像区域;(2)基于图像观测模型,对眼底图像进行归一化增强;(3)采用非血管结构抑制算子,同时结合滞后多阈值处理技术,实现眼底图像血管提取与分割;(4)设定2倍主血管宽度的垂直窗口,在血管分割图上沿水平方向滑动并计算各水平位置x的血管散布度D(x),从而获得水平投影的散布度曲线,该曲线的最小值点确定为视盘水平坐标x
1.一种综合血管分布和视盘外观特性的视盘投影定位方法,包括以下步骤:(1)取原彩色眼底图像中红色通道分量最亮的像素强度值的5%作为阈值,并根据此阈值分割出前景区域的二值图像,对二值图像进行形态学腐蚀操作,得到掩膜模板,通过掩膜处理得到眼底图像的感兴趣区域;(2)根据图像观察模型,对掩膜处理后的感兴趣区域进行照度和对比度的归一化处理,得到增强的视网膜眼底图像;(3)对步骤(2)得到的增强视网膜眼底图像采用Gabor滤波非血管结构抑制算子,同时结合滞后多阈值处理技术,实现眼底图像血管提取与分割,得到血管分割图;(4)设置一个垂直窗口,窗口的高度为图像高度,宽度为2倍的主血管宽度,沿血管分割图的水平方向从左到右滑动;在窗口内计算每一个水平位置x处的血管散布度值D(x),根据血管散布度值绘出水平投影曲线,找到该一维曲线的最小值点确定为视盘的水平坐标x_od;(5)设置一个矩形窗口,窗口的长度和宽度为视盘直径的大小,且其中心的水平坐标为已确定的视盘水平坐标x_od,分别在原灰度图和初相位ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="171"/>的Gabor滤波响应图上沿垂直方向从上到下滑动窗口,并估计窗口内平均亮度和平均滤波响应值,将平均亮度值和平均滤波响应值相乘作为相应垂直位置y处的纵向投影值,根据纵向投影值绘制投影曲线,找到最大值点的坐标为视盘的垂直坐标。
2.根据权利要求1所述的综合血管分布和视盘外观特性的视盘投影定位方法,所述步骤(1)中的形态学腐蚀操作是采用9×9方形结构元素对二值图像进行腐蚀运算。
3.根据权利要求1所述的综合血管分布和视盘外观特性的视盘投影定位方法,所述步骤(2)的具体步骤为:1将图像分为大小为s的分块S_i,s为1/6~1/8的M×N,M×N为原始图像尺寸,对每个分块S_i计算该分块中的灰度均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="136"/>和标准差ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="168"/>再采用双三次插值方法求得整幅图像中每点的均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="58"/>与标准差ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="84"/>2通过Manhabolios距离d_M(x,y),判断某像素(x,y)是否属于背景区域,如果其距离邻域均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="55"/>小于某阈值t,阈值t=1,即 d M ( x , y ) = | I ( x , y ) - μ ^ N ( x , y ) σ ^ N ( x , y ) | < t - - - ( 11 ) ]]> 则该像素属于背景区域,否则属于前景区域;3在分割出来的背景区域对照度和对比度进行估计,每一个分块中采样点(s_x,s_y)的照度ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="156"/>和对比度ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="175"/>由分块中背景像素点的均值和标准差来估计获得,其余非采样点和非背景像素点的照度和对比度漂移因子则通过双三次插值方法获得;4根据式(7),实现图像的归一化增强, I 0 ^ ( x , y ) = I ( x , y ) - L ^ ( x , y ) C ^ ( x , y ) - - - ( 7 ) ]]> I(x,y)为原始观测图像在(x,y)像素点的灰度值,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="129"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="133"/>分别为该点估算的照度和对比度漂移因子,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="152"/>是归一化后的灰度值。
4.根据权利要求1所述的综合血管分布和视盘外观特性的视盘投影定位方法,所述步骤(3)具体步骤为:1选用12个不同方向上的Gabor滤波模板分别对图像进行滤波,每个相邻方向相差π/12,滤波后的图像用ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="257"/>表示: ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> 其中,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="240"/>为Gabor函数,f代表输入图像,保留每个像素在12个不同方向的最大响应, ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> H ( z ) = 0 , z < 0 z , z ≥ 0 - - - ( 15 ) ]]> 2针对Gabor函数在不同相位ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="35"/>时,对不同特征模式的检测特性,采用式(16)对明条纹及非血管边缘进行抑制处理, v λ , δ ( x , y ) = H ( r λ , δ , π AH ( x , y ) - α 1 r λ , δ , 0 AH ( x , y ) - α 2 r λ , δ , π / 2 AH ( x , y ) ) - - - ( 16 ) ]]> 其中ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="700"/>分别表示ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="135"/>0,π/2时的滤波响应,α₁,α₂分别表示对明条纹及非血管边缘的抑制分量ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="464"/>的抑制程度;3选取多组阈值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="700"/>考虑到血管的连通性,相邻两组阈值区间应该重叠,即ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="469"/>选定阈值对ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="60"/>对滤波抑制处理后图像进行单对阈值二值化分割,得到相应的二值化图像ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="208"/> b λ , δ i ( x , y ) = H T i ( v λ , δ ( x , y ) ) , i ∈ { 1 , . . . , N i } - - - ( 17 ) ]]> HT表示标准滞后阈值处理;4将每一个二值图像ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="180"/>表示成若干连通分量的并集, b λ , δ i = ∪ k C λ , δ i , k - - - ( 18 ) ]]> 其中,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="314"/>表示阈值对选取为ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="29"/>时的二值化图像中第k个连通分量。利用形态学膨胀操作,选用2×2方形结构元素q₂,对ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="154"/>的每个分支ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="78"/>进行膨胀处理,得到: D λ , δ i , k = C λ , δ i , k ⊕ q 2 - - - ( 19 ) ]]> 5最终的血管结构V_λ,δ可以通过下述公式依序处理获得, B λ , δ ( N t , N t - 1 ) = ∪ ( b λ , δ N t ⊕ q 2 ) ∩ D λ , δ N t - 1 , k ≠ φ ( b λ , δ N t ∪ C λ , δ N t - 1 , k ) - - - ( 20 ) ]]> B λ , δ ( N t - 1 , N t - 2 ) = ∪ ( B λ , δ ( N t , N t - 1 ) ⊕ q 2 ) ∩ D λ , δ N t - 2 , k ≠ φ ( B λ , δ N t , N t - 1 ∪ C λ , δ N t - 2 , k ) - - - ( 21 ) ]]> ... V λ , δ = B λ , δ ( 2,1 ) = ∪ ( B λ , δ ( 3,2 ) ⊕ q 2 ) ∩ D λ , δ 1 , k ≠ φ ( B λ , δ ( 3,2 ) ∪ C λ , δ 1 , k ) - - - ( 22 ) ]]> 6由于血管在图像中呈现连通的网络结构,考虑去除所有长度小于10个像素的分支。
5.根据权利要求1所述的综合血管分布和视盘外观特性的视盘投影定位方法,所述的步骤(4)中血管散布度值D(x)的计算公式为: D ( x ) = ( - 1 ) Σ i = 1 n x p i log 2 p i max i { p i } - - - ( 23 ) ]]> 上述公式中ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="176"/>为该位置垂直窗内第i段血管连通分量所占比重,m_i表示第i段连通血管像素数,M表示水平滑动窗口内血管总的像素个数,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="195"/>表示最大连通血管段的像素个数,n_x为该位置水平滑动窗口内总的血管连通数目。
6.根据权利要求1所述的综合血管分布和视盘外观特性的视盘投影定位方法,所述步骤(5)的具体步骤为:1定义窗口矩形窗口W_v,W_v的长度和宽度为视盘直径大小,且其中心的水平坐标已确定的视盘水平坐标x_od处,纵坐标y随窗口上、下滑动而变化;2沿垂直方向从上到下滑动窗口W_v,并分别在原图和Gabor滤波在θ=90°相位的响应图上估计窗口内平均亮度和平均滤波响应值,根据式(24),将平均亮度值和平均滤波响应值相乘作为相应垂直位置y处的纵向投影值,f(y)=I_N(x_od,y)*g_N(x_od,y)xa0xa0xa0xa0xa0xa0xa0xa0xa0(24)I_N(x_od,y)为以(x_od,y)为中心的矩形窗邻域的平均亮度值,g_N(x_od,y)为该邻域内平均Gabor滤波响应估计值;3根据f(y)绘制纵向投影曲线,找到最大值点的坐标为视盘的垂直坐标y_od。技术领域
本发明涉及一种视网膜眼底图像中的视盘进行自动定位的方法,特别涉及一种综合血管分布和视盘外观特性的视盘投影定位方法。
背景技术
视盘是视网膜的主要生理结构之一,视神经和血管从该区域进入眼部并向周边辐射延伸,在视网膜眼底图像中表现为圆形亮黄色区域,同时包含大量较粗的血管。研究人员对视盘的自动定位一直非常关注,因为视盘的准确定位有利于眼底病变的诊断,例如:很多视盘分割算法需要预先给出初始种子点;由于视盘和黄斑(视觉最敏锐区域)之间的距离近似为常数,当视盘位置确定后,确定黄斑区域也就容易了;视盘和大尺寸的亮黄色病变(溢出物病变)容易混扰,视盘的定位可以将其从候选病变对象中予以排除。目前视盘定位方法大致分为以下几类:1、基于视盘外观特性的方法。该类方法通常采用视盘的外观如亮度、对比度、形状信息作为定位特征,一般将亮度最大或对比度最强的圆形区域的中心作为视盘的参考位置。单纯利用外观特性,在质量较好的正常眼底图像中的视盘定位的成功率较高,但是在病变图像中,由于病变造成视盘外观改变,或其它大面积亮色病变区域干扰,视盘定位容易出错。2、基于血管特性的方法。这主要是根据视盘是血管进入眼部的起始区域和血管分支汇合点,且该区域血管最粗,血管密度最大。基于血管特性检测的视盘定位算法一般都比较复杂,而且耗时,同时在低质量或病变图像中,血管特征检测仍是一个比较困难的问题。3、基于视盘区域的多种特征的方法。该类方法综合利用视盘的外观特性和血管特性,定位精度相对较高。但是该类方法需要建立较复杂模型,其算法复杂,不便于实时应用。近年来,国内外众多学者针对视盘定位进行了大量的研究,尽管如此,绝大部分研究在视盘定位的准确率和速度两方面难以兼顾,导致视盘定位方法很难适应实际系统的实时处理要求。
发明内容
为了解决现有视盘定位存在的上述技术问题,本发明提供了一种综合血管分布和外观特性的视盘投影定位方法。本发明综合利用了视盘区域的血管分布和外观特性信息,同时由于将复杂的二维定位问题转换为相对简单的一维定位问题,算法实现简单,成功率高,具有良好的鲁棒性,综合准确率和处理速度两者来看,本发明符合实际系统实时处理要求。本发明解决上述技术问题的技术方案包括以下步骤:(1)取原彩色眼底图像中红色通道分量最亮的像素强度值的5%作为阈值,并根据此阈值分割出前景区域的二值图像,对二值图像进行形态学腐蚀操作,得到掩膜模板,通过掩膜处理得到眼底图像的感兴趣区域;(2)根据图像观察模型,对掩膜处理后的感兴趣区域进行照度和对比度的归一化处理,得到增强的视网膜眼底图像;(3)对步骤(2)得到的增强视网膜眼底图像采用Gabor滤波非血管结构抑制算子,同时结合滞后多阈值处理技术,实现眼底图像血管提取与分割,得到血管分割图;(4)设置一个垂直窗口,窗口的高度为图像高度,宽度为2倍的主血管宽度,沿血管分割图的水平方向从左到右滑动;在窗口内计算每一个水平位置x处的血管散布度值D(x),根据血管散布度值绘出水平投影曲线,找到该一维曲线的最小值点确定为视盘的水平坐标x_od;(5)设置一个矩形窗口,窗口的长度和宽度为视盘直径的大小,且其中心的水平坐标为已确定的视盘水平坐标x_od,分别在原灰度图和初相位ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="171"/>的Gabor滤波响应图上沿垂直方向从上到下滑动窗口,并估计窗口内平均亮度和平均滤波响应值,将平均亮度值和平均滤波响应值相乘作为相应垂直位置y处的纵向投影值,根据纵向投影值绘制投影曲线,找到最大值点的坐标为视盘的垂直坐标。上述的综合血管分布和视盘外观特性的视盘投影定位方法中,所述步骤(1)中的形态学腐蚀操作是采用9×9方形结构元素对二值图像进行腐蚀运算。上述的综合血管分布和视盘外观特性的视盘投影定位方法中,所述步骤(2)的具体步骤为:1将图像分为大小为s的分块S_i,s为1/6~1/8的M×N,M×N为原始图像尺寸,对每个分块S_i计算该分块中的灰度均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="136"/>和标准差ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="168"/>再采用双三次插值方法求得整幅图像中每点的均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="58"/>与标准差ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="84"/>2通过Manhabolios距离d_M(x,y),判断某像素(x,y)是否属于背景区域,如果其距离邻域均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="55"/>小于某阈值t,阈值t=1,即 d M ( x , y ) = | I ( x , y ) - μ ^ N ( x , y ) σ ^ N ( x , y ) | < t - - - ( 11 ) ]]> 则该像素属于背景区域,否则属于前景区域;3在分割出来的背景区域对照度和对比度进行估计,每一个分块中采样点(s_x,s_y)的照度ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="156"/>和对比度ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="175"/>由分块中背景像素点的均值和标准差来估计获得,其余非采样点和非背景像素点的照度和对比度漂移因子则通过双三次插值方法获得;4根据式(7),实现图像的归一化增强, I 0 ^ ( x , y ) = I ( x , y ) - L ^ ( x , y ) C ^ ( x , y ) - - - ( 7 ) ]]> I(x,y)为原始观测图像在(x,y)像素点的灰度值,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="129"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="133"/>分别为该点估算的照度和对比度漂移因子,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="151"/>是归一化后的灰度值。上述的综合血管分布和视盘外观特性的视盘投影定位方法中,所述步骤(3)具体步骤为:1选用12个不同方向上的Gabor滤波模板分别对图像进行滤波,每个相邻方向相差π/12,滤波后的图像用ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="257"/>表示: 其中,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="240"/>为Gabor函数,f代表输入图像,保留每个像素在12个不同方向的最大响应, ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> H ( z ) = 0 , z < 0 z , z ≥ 0 - - - ( 15 ) ]]> 2针对Gabor函数在不同相位ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="35"/>时,对不同特征模式的检测特性,采用式(16)对明条纹及非血管边缘进行抑制处理, v λ , δ ( x , y ) = H ( r λ , δ , π AH ( x , y ) - α 1 r λ , δ , 0 AH ( x , y ) - α 2 r λ , δ , π / 2 AH ( x , y ) ) - - - ( 16 ) ]]> 其中ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="700"/>分别表示ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="135"/>0,π/2时的滤波响应,α₁,α₂分别表示对明条纹及非血管边缘的抑制分量ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="464"/>的抑制程度;3选取多组阈值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="700"/>考虑到血管的连通性,相邻两组阈值区间应该重叠,即ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="469"/>选定阈值对ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="60"/>对滤波抑制处理后图像进行单对阈值二值化分割,得到相应的二值化图像ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="207"/> b λ , δ i ( x , y ) = H T i ( v λ , δ ( x , y ) ) , i ∈ { 1 , . . . , N i } - - - ( 17 ) ]]> HT表示标准滞后阈值处理;4将每一个二值图像ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="180"/>表示成若干连通分量的并集, b λ , δ i = ∪ k C λ , δ i , k - - - ( 18 ) ]]> 其中,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="314"/>表示阈值对选取为ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="29"/>时的二值化图像中第k个连通分量。利用形态学膨胀操作,选用2×2方形结构元素q₂,对ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="153"/>的每个分支ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="78"/>进行膨胀处理,得到: D λ , δ i , k = C λ , δ i , k ⊕ q 2 - - - ( 19 ) ]]> 5最终的血管结构V_λ,δ可以通过下述公式依序处理获得, B λ , δ ( N t , N t - 1 ) = ∪ ( b λ , δ N t ⊕ q 2 ) ∩ D λ , δ N t - 1 , k ≠ φ ( b λ , δ N t ∪ C λ , δ N t - 1 , k ) - - - ( 20 ) ]]> B λ , δ ( N t - 1 , N t - 2 ) = ∪ ( B λ , δ ( N t , N t - 1 ) ⊕ q 2 ) ∩ D λ , δ N t - 2 , k ≠ φ ( B λ , δ N t , N t - 1 ∪ C λ , δ N t - 2 , k ) - - - ( 21 ) ]]> ... V λ , δ = B λ , δ ( 2,1 ) = ∪ ( B λ , δ ( 3,2 ) ⊕ q 2 ) ∩ D λ , δ 1 , k ≠ φ ( B λ , δ ( 3,2 ) ∪ C λ , δ 1 , k ) - - - ( 22 ) ]]> 6由于血管在图像中呈现连通的网络结构,考虑去除所有长度小于10个像素的分支。上述的综合血管分布和视盘外观特性的视盘投影定位方法,所述的步骤(4)中血管散布度值D(x)的计算公式为: D ( x ) = ( - 1 ) Σ i = 1 n x p i log 2 p i max i { p i } - - - ( 23 ) ]]> 上述公式中ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="176"/>为该位置垂直窗内第i段血管连通分量所占比重,m_i表示第i段连通血管像素数,M表示水平滑动窗口内血管总的像素个数,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="195"/>表示最大连通血管段的像素个数,n_x为该位置水平滑动窗口内总的血管连通数目。上述的综合血管分布和视盘外观特性的视盘投影定位方法中,所述步骤(5)的具体步骤为:1定义窗口矩形窗口W_v,W_v的长度和宽度为视盘直径大小,且其中心的水平坐标已确定的视盘水平坐标x_od处,纵坐标y随窗口上、下滑动而变化;2沿垂直方向从上到下滑动窗口W_v,并分别在原图和Gabor滤波在θ=90°相位的响应图上估计窗口内平均亮度和平均滤波响应值,根据式(24),将平均亮度值和平均滤波响应值相乘作为相应垂直位置y处的纵向投影值,f(y)=I_N(x_od,y)*g_N(x_od,y)xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(24)I_N(x_od,y)为以(x_od,y)为中心的矩形窗邻域的平均亮度值,g_N(x_od,y)为该邻域内平均Gabor滤波响应估计值;3根据f(y)绘制纵向投影曲线,找到最大值点的坐标为视盘的垂直坐标y_od。与现有技术相比,本发明的技术效果在于:1、现有血管提取方法很少考虑到视网膜病变后,由于受病变区域边界干扰,很难正确提取血管结构。视盘、出血或其他病变区域的边界具有类似血管的特性,这些边界在病变视网膜眼底图像中大量存在。如何在血管特征提取阶段抑制非血管区域边缘效应是本发明着重解决的问题之一。本发明利用不同相位下,Gabor滤波函数对不同特征模式(暗条纹、亮条纹和边缘)的响应特性,提出一种基于Gabor滤波的非血管结构抑制算子,同时结合滞后多阈值处理技术,最终实现眼底图像血管结构提取。本发明所用血管提取方法能够很大程度上抑制病变图像中的非血管像素,并较完整的保留了血管结构,极大提高了血管提取的准确性。2、本发明综合利用了血管分布和视盘外观特性。首先根据血管的分布特性(视盘区域血管沿垂直方向伸展,联通性好,分布集中)确定视盘的水平坐标,进而利用视盘外观特性(亮度、边缘梯度)定位垂直坐标,由于将复杂的二维定位问题转换为相对简单的两个一维定位问题,算法实现简单,成功率高,具有良好的鲁棒性。
附图说明
图1是本发明的流程图。图2是本发明中彩色眼底图像掩膜模板示意图。图3是本发明中的眼底图像归一化示意图。图4是不同相位是Gabor函数空间模板示意图。图5是某眼底彩色图像及其和对应的血管分布图。图6是在眼底灰度图像及其ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="171"/>相位Gabor滤波图像上的矩形窗示意图。图7是某眼底图像的坐标投影定位示意图。
具体实施方式
以下将结合附图和实施例对本发明做进一步详细说明。如图1所示,本发明的视盘投影定位方法,其具体流程为:1、眼底图像掩膜处理。如图2(a)所示,视网膜眼底图像通常包含暗色背景和视网膜眼底成像区域,其中视网膜眼底成像区域是眼科专家的感兴趣区域,为了防止背景和边界的干扰,有必要获得感兴趣区域的掩膜模板。由于彩色眼底图的红色分量接近饱和,最能反映成像时照明情况,因此取原彩色眼底图像中红色通道分量I_R最亮的像素强度值的5%作为分割阈值t_b,并根据此阈值分割前景区域的二值图像,t_b=0.05*max(I_R)xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(1) I bin ( x , y ) = 1 , I R ( x , y ) > t b 0 , I R ( x , y ) ≤ t b - - - ( 2 ) ]]> 为尽可能排除边界的干扰,有必要对二值图I_bin进行腐蚀操作,本发明采用9×9方形结构q₉元素对I_bin进行腐蚀运算,获得最终的掩膜模板I_mask,参见图2(b)。 ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="46"/>表示腐蚀运算。采用掩膜操作不难得到图像的感兴趣区域。2、眼底图像归一化增强。由于图像获取过程中照明、曲面反射、视角、眼睛运动的影响,眼底图像普遍存在非均匀照度和对比度现象,因此眼底图像的归一化是后续处理的必要环节。原始视网膜眼底图像可看作理想背景图像和前景图像的叠加组合,在图像观测(获取)过程中,图像照度和对比度会在原始图像基础上产生形变,这种变形可以用如下的观测模型描述:I(x,y)=C(x,y)*I⁰(x,y)+L(x,y)xa0xa0xa0xa0xa0xa0xa0(4) I 0 ( x , y ) = I b 0 ( x , y ) + I f 0 ( x , y ) - - - ( 5 ) ]]> I⁰是原始眼底图像,是不包含非均匀照度和对比度影响的理想眼底图像。ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="39"/>是原始背景图像,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="42"/>是原始前景图像。L(x,y)和C(x,y)分别是获取过程中像素点(x,y)的照度漂移因子和对比度漂移因子,I是最终观测到的图像。根据式(4),可以获得原始眼底图像, I 0 ( x , y ) = I ( x , y ) - L ( x , y ) C ( x , y ) - - - ( 6 ) ]]> 真实的照度漂移因子和对比度因子通常是未知的,只能从观测图像中对其值进行估计。因此式(6)变成 I 0 ^ ( x , y ) = I ( x , y ) - L ^ ( x , y ) C ^ ( x , y ) - - - ( 7 ) ]]> ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="40"/>是对原始眼底图像I⁰的估计。为获得ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="68"/>必须估计出照度漂移L和对比度漂移C。由于前景部分特性差异很大,因此前景部分ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="42"/>的特性难以表达,另一方面,背景ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="40"/>变化平缓,可以通过正态分布规律来描述, I b 0 ( x , y ) ~ N ( μ b , σ b ) - - - ( 8 ) ]]> 其中μ_b表示背景的理想均匀照度值,标准差σ_b则反映了视网膜在空间域的自然特性变化。因为照度的不规则并不会造成图像的剧烈变化,意味着照度漂移L和对比度漂移C都主要集中在频谱的低频成分,所以背景图像是估计L和C的主要依据。为实现漂移因子的估计,不妨将原始图像分割为前景图像和背景图像,这样可以将注意力放在背景部分,在只考虑背景区域时,可以令相关前景区域ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="155"/>因此式(4)可以简化为: I ( x , y ) = C ( x , y ) * I 0 ( x , y ) + L ( x , y ) ]]> = C ( x , y ) * [ I b 0 ( x , y ) + I f 0 ( x , y ) ] + L ( x , y ) - - - ( 9 ) ]]> = C ( x , y ) * I b 0 ( x , y ) + C ( x , y ) * I f 0 ( x , y ) + L ( x , y ) ]]> ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="37"/>是背景区域像素构成的集合。显而易见,背景符合正态分布特性: I b 0 ( x , y ) ~ N ( μ b , σ b ) ~ N ( L ( x , y ) , C ( x , y ) ) - - - ( 10 ) ]]> 从上面模型的分析可知,可以通过观测图像的背景区域ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="434"/>的均值和标准差来估计背景中的ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="114"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="146"/>并最终根据式(7)恢复原始图像。整个过程可以分为两步:第一步首先是背景提取,然后第二步通过背景区域估计照度和对比度漂移因子ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="112"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="149"/>由于RGB图像中的绿色通道对比度最强,相对其它两通道能更好的反映眼底的特性,在此后续处理都是在绿色通道上进行的。其具体实现过程如下:1将图像分为大小为s(约为1/6~1/8的M×N,M×N为原始图像尺寸)的分块S_i,然后对每个分块S_i计算该分块中的均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="136"/>和标准差ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="168"/>最后采用双三次插值方法求得整幅图像中每点的ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="58"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="59"/>值;2通过Manhabolios距离d_M(x,y),判断某像素(x,y)是否属于背景区域,如果其距离邻域均值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="55"/>小于某阈值t(阈值t=1),即 d M ( x , y ) = | I ( x , y ) - μ ^ N ( x , y ) σ ^ N ( x , y ) | < t - - - ( 11 ) ]]> 则该像素属于背景区域,否则属于前景区域;3在分割出来的背景区域对照度和对比度进行估计。由于背景只是整幅眼底图像的一部分,为避免逐点计算,采用类似第1步的样本点取样估计。每一个分块中采样点(s_x,s_y)的照度ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="157"/>和对比度ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="170"/>由分块中背景像素点的均值和标准差来估计获得,其余非采样点和非背景像素点的照度和对比度漂移因子则通过双三次插值方法获得。最终,获得整幅眼底图像每一个像素点的ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="137"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="172"/>4根据式(7),实现图像的归一化增强。图3(a-e)是某示例图的归一化处理过程示意图,其中图3(a)为原始图像,图3(b)为图像的背景估计图,图3(c)和图3(d)分别是亮度漂移因子和对比度漂移因子的估计示意图,图3(e)是最终的归一化彩图。3、血管提取与分割。虽然已有很多文献提出了有效的视网膜血管提取方法,但是很少文献提到病变图像中的血管提取问题,病变区域和视盘的边缘在具有强对比度时,其匹配滤波响应(例如Gabor滤波器)比大多数弱、细血管更强,因此在采用简单阈值分割时,相关边缘常常被误判为血管。为获得尽可能准确完整的血管结构,本发明提出一种非血管抑制算子,同时结合多滞后阈值技术实现血管的提取。具体流程如下:1由于血管为狭长暗色条状或线状结构,其横截面近似于高斯分布,我们采用具有可调整视场和方向的Gabor滤波器实现血管特征的提取, ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> x ~ = x cos θ + y sin θ ]]> y ~ = - x sin θ + y cos θ ]]> 其中,γ是空间长宽比,决定了接受视场的椭圆度,本发明中设为常数0.5,标准差δ确定接受视场的尺寸。参数λ是余弦因子的波长。比值δ/λ决定空间频域的带宽,影响到视场中平行的兴奋性的和抑制性的条纹区域数量,文中取为定值δ/λ=0.56。角度参数θ,θ∈[0,π),为预定的滤波方向。相位偏移参数ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="59"/>ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="220"/>确定ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="240"/>函数的对称性。Gabor滤波器的响应ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="226"/>由Gabor函数ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="240"/>和输入图像f的卷积计算获得。选用12个不同方向上({θ_i=0,π/12,2π/12,...11π/12}|,每个相邻方向相差π/12)的Gabor滤波模板分别对图像进行滤波,滤波后的图像用ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="256"/>表示: ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> 保留每个像素在12个不同方向的最大响应, ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> H ( z ) = 0 , z < 0 z , z ≥ 0 - - - ( 15 ) ]]> 其中,H为半波整定函数。在理想条件下,血管中心像素在相位偏移参数ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="115"/>时具有最大响应,在ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="104"/>时响应接近于0。而边缘在ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="189"/>或π/2时有最大响应值。2血管在视网膜灰度图像中呈现暗条纹结构。当Gabor滤波函数ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="109"/>时,核函数空间结构如图4(a)所示,此时Gabor滤波函数对暗条纹有最大响应;当ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="100"/>时,核函数空间结构如图4(b),此时Gabor滤波函数对亮条纹响应最大;当ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="169"/>时,核函数空间结构如图4(c),此时边缘有最大响应。针对Gabor函数在不同相位ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="35"/>时,对不同特征模式的检测特性,采用式(16)对明条纹及非血管边缘进行抑制处理, v λ , δ ( x , y ) = H ( r λ , δ , π AH ( x , y ) - α 1 r λ , δ , 0 AH ( x , y ) - α 2 r λ , δ , π / 2 AH ( x , y ) ) - - - ( 16 ) ]]> 其中ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="700"/>分别表示ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="109"/>(暗条纹),0(亮条纹),π/2(边缘)时的滤波响应,α₁,α₂分别表示对明条纹及非血管边缘的抑制分量ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="465"/>的抑制程度。在实际提取过程中,选择α₁=0.5、α₂=1,此时抑制算法对亮条纹以及非血管边缘等区域均有较好的综合抑制效果。3通过简单阈值分割可以提取主要的血管结构。尽管大部分非血管像素在抑制算子处理后得到有效抑制,但仍存在一些散布的点状和短线状残留像素,为去除这些残留信息,在标准滞后阈值方法基础上,本发明提出了一种多滞后阈值血管分割方法。标准滞后阈值包括一对阈值t_l和t_h,t_lh。Gabor响应值大于t_h的像素在输出二值血管图时指定为1,而Gabor响应值小于t_l的像素被指定为0,响应值在t_l和t_h之间的像素,如果能够找到一条联接通道可以和某一响应值大于t_h的点连通,通道上的像素响应值大于t_l,则该像素也被指定为1。基于标准滞后阈值技术,选取多组阈值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="700"/>考虑到血管的连通性,相邻两组阈值区间应该重叠,即ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="474"/>例如,在实现例中采用4对阈值ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="74"/>i=1,...,4,阈值参数对选取为:[0.02,0.04]*max(v_λ,δ),[0.03,0.05]*max(v_λ,δ),[0.04,0.06]*max(v_λ,δ),[0.05,0.08]*max(v_λ,δ),max(v_λ,δ)为整幅图像中v_λ,δ的最大值。选定阈值对ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="62"/>对滤波抑制处理后图像进行单对阈值二值化分割,得到相应的二值化图像ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="206"/> b λ , δ i ( x , y ) = H T i ( v λ , δ ( x , y ) ) , i ∈ { 1 , . . . , N i } - - - ( 17 ) ]]> HT表示标准滞后阈值处理。4将每一个二值图ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="180"/>表示成若干连通分量的并集, b λ , δ i = ∪ k C λ , δ i , k - - - ( 18 ) ]]> 其中,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="314"/>表示阈值对选取为ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="25"/>和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="29"/>时的二值化图像中第k个连通分量。利用形态学膨胀操作,选用2×2方形结构元素q₂,对ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="154"/>的每个分支ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="78"/>进行膨胀处理,得到: D λ , δ i , k = C λ , δ i , k ⊕ q 2 - - - ( 19 ) ]]> 5最终的血管结构V_λ,δ可以通过下述公式依序处理获得, B λ , δ ( N t , N t - 1 ) = ∪ ( b λ , δ N t ⊕ q 2 ) ∩ D λ , δ N t - 1 , k ≠ φ ( b λ , δ N t ∪ C λ , δ N t - 1 , k ) - - - ( 20 ) ]]> B λ , δ ( N t - 1 , N t - 2 ) = ∪ ( B λ , δ ( N t , N t - 1 ) ⊕ q 2 ) ∩ D λ , δ N t - 2 , k ≠ φ ( B λ , δ N t , N t - 1 ∪ C λ , δ N t - 2 , k ) - - - ( 21 ) ]]> ..x7f. V λ , δ = B λ , δ ( 2,1 ) = ∪ ( B λ , δ ( 3,2 ) ⊕ q 2 ) ∩ D λ , δ 1 , k ≠ φ ( B λ , δ ( 3,2 ) ∪ C λ , δ 1 , k ) - - - ( 22 ) ]]> 6由于血管在图像中呈现连通的网络结构,考虑去除所有长度小于10个像素的分支。本发明采用的血管提取方法提取的血管主干结构完整、干净,能够排除大部分非血管像素。4、基于血管分布特性的水平坐标定位。图5(b)是某眼底图像(图5(a))的人工标记的血管二值图,在限定的宽度约为2倍的主血管宽度的垂直窗内,不难观察到在不同水平位置,血管分布有显著差异。在视盘位置(图5(b)),血管联通数目少,而且分布非常集中,局部血管密度大,而在其它位置,血管连通数目多而且分布相对零散。根据这一分布特性,视盘水平坐标定位的具体流程如下:1定义窗口W_h(高度为图像高度,宽度大约为2倍的主血管宽度),且其中心定位在水平坐标点x处;由于血管的连通性,血管主干结构对于血管散布度的计算更为重要,为尽可能排除弱、细血管干扰,计算血管散布度前可以先将定义的垂直窗口内一些较短的小血管段去掉,根据经验,阈值th=30,即忽略掉连通血管像素小于30的小段血管。2在血管分割图的水平方向从左到右滑动窗口W_h,并计算每个窗口内的血管散布度值D(x),作为水平位置x处的投影值; D ( x ) = ( - 1 ) Σ i = 1 n x p i log 2 p i max i { p i } - - - ( 23 ) ]]> 其中D(x)表示水平位置x处的血管散布度,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="176"/>为该位置垂直窗内第i段血管连通分量所占比重,m_i表示第i段连通血管像素数,M表示垂直窗口血管总的像素个数。ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="196"/>表示最大连通血管段的像素个数,可以一定程度衡量血管的局部密度,n_x为该位置垂直窗内总的血管连通数目。一般情况下,包含视盘的垂直窗口区域对应的散布度要小于其它水平位置的散布度。垂直窗沿水平方向滑动,可以得到各水平位置的血管散布度值,从而绘出其在水平方向的一维投影曲线如图7(c),无血管或血管像素占垂直窗口面积少于0.01%的位置的血管散布度预设为无穷大值。3水平方向的一维投影曲线的最小值点即可确定为视盘的水平坐标x_od;5、基于视盘外观特性的垂直坐标定位。由于视盘区域的亮度和边缘梯度值相对要大于其他区域,在精准定位视盘的水平坐标后,可利用这一外观特性进一步定位视盘垂直坐标。具体做法为:1选择视盘大小的矩形窗W_v(长度和宽度为视盘直径大小),且其中心的水平坐标已确定的视盘水平坐标x_od处,纵坐标y随窗口上、下滑动而变化;2矩形窗W_v分别在原灰度图(图6(a))和ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="170"/>相位Gabor滤波响应图(图6(b))上沿垂直方向上下滑动。并分别在原图和Gabor滤波在ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="171"/>相位的响应图上估计窗口内平均亮度和平均滤波响应值,根据式(24),将平均亮度值和平均滤波响应值相乘作为相应垂直位置y处的纵向投影特征值f(y),f(y)=I_N(x_od,y)*g_N(x_od,y)xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(24)其中I_N(x_od,y)为以(x_od,y)为中心的矩形窗邻域的平均亮度值,g_N(x_od,y)为该邻域内Gabor滤波在ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="171"/>相位时的平均滤波响应估计值,反映了领域内的边缘梯度信息。由于视盘区域的亮度和血管边缘信息都要大于其它区域,因此f(y)值大的位置为视盘的可能性大,同样在矩形窗上下滑动时,可以投影得到相应位置的反映f(y)变化情况的一维曲线,如图7(b)。3在一维f(y)投影曲线上找到最大值点的坐标为视盘的垂直坐标y_od。以上所述仅是本发明的优选实施方式,本发明的保护范围并不仅局限于上述实施例,凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理前提下的若干改进和润饰,应视为本发明的保护范围。